home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / oper_sys / weyl / weyl_lsp.lha / grobner.lisp < prev    next >
Encoding:
Text File  |  1991-11-13  |  10.3 KB  |  311 lines
  1. ;;; -*- Mode:Lisp; Package:Weyli; Base:10; Lowercase:T; Syntax:Common-Lisp -*-
  2. ;;; ===========================================================================
  3. ;;;                  Expanded Polynomials
  4. ;;; ===========================================================================
  5. ;;; (c) Copyright 1989, 1991 Cornell University
  6.  
  7. ;;; $Id: grobner.lisp,v 2.6 1991/11/13 16:15:11 rz Exp $
  8.  
  9. (in-package "WEYLI")
  10.  
  11. (defmacro with-grobner-operations (grobner-basis &body body)
  12.   `(with-slots (compare-function ring generators undones reducibles possibles)
  13.          ,grobner-basis
  14.      (macrolet ((e> (a b) `(%funcall compare-function ,a ,b))
  15.         (e< (a b) `(%funcall compare-function ,b ,a)))
  16.        ,@body)))
  17.  
  18. ;; Grobner calculations are done within the context of an instance of
  19. ;; the Grobner-Basis flavor.  Each instance has its own variable list
  20. ;; and flags sets.  At any time the user can add polynomials or
  21. ;; extract information from the structure.
  22.  
  23. (defclass grobner-basis ()
  24.   ((ring :initarg :ring)
  25.  
  26.    ;; The exponent comparison function
  27.    (compare-function :initform #'lexical->
  28.              :initarg :compare-function
  29.              :reader compare-function)
  30.  
  31.    ;; the current list of generators
  32.    (generators :initform ())
  33.  
  34.    (undones :initform ())
  35.    ;; A list of triples pairs (lt f g), lt(f)<=lt(g), of elements of
  36.    ;; GENERATORS such that if any pair is not in the list, its s-poly
  37.    ;; is guaranteed to be writable as a linear combination of
  38.    ;; GENERATORS, with smaller s-degs
  39.    
  40.    (reducibles :initform nil)
  41.    (possibles :initform nil)    
  42.    ))
  43.  
  44. (defmethod print-object ((grob-struct grobner-basis) stream)
  45.   (let ((gens (relations grob-struct)))
  46.     (format stream "(~S~{, ~S~})"
  47.         (first gens) (rest gens))))
  48.  
  49. (defun check-same-domain (exprs)
  50.   (let ((domain (domain-of (first exprs))))
  51.     (loop for exp in (rest exprs)
  52.       do (unless (eql domain (domain-of exp))
  53.            (return nil))
  54.       finally (return domain))))
  55.  
  56. (defun make-ideal (&rest polys)
  57.   (let ((domain (check-same-domain polys))
  58.     ideal)
  59.     (when (null domain)
  60.       (error "Not all polynomials are from the same domain:~%~S"
  61.          polys))
  62.     (unless (typep domain 'polynomial-ring)
  63.       (error "~S is not a polynomial ring"))
  64.     (cond ((typep (coefficient-domain domain) 'field)
  65.        (setq ideal (make-instance 'grobner-basis :ring domain)))
  66.       (t (error "Can't deal with polynomials not over fields: ~S"
  67.             domain)))
  68.     (loop for p in polys
  69.       do (add-relation ideal p))
  70.     ideal))
  71.  
  72. (defmethod (setf compare-function) (new-function (grob grobner-basis))
  73.   (with-slots (ring compare-function generators reducibles possibles) grob
  74.     (unless (eql compare-function new-function) 
  75.       (flet ((convert-list (list)
  76.            (loop for poly in list
  77.              collect
  78.              (sort poly
  79.                #'(lambda (a b)
  80.                    (%funcall new-function
  81.                      (first a) (first b)))))))
  82.     (setq generators (convert-list generators))
  83.     (setq reducibles (convert-list reducibles))
  84.     (setq possibles (convert-list possibles))
  85.     (setq compare-function new-function)))
  86.     grob))    
  87.  
  88. (defmethod add-relation ((grob-struct grobner-basis) (relation mpolynomial))
  89.   (with-slots (ring compare-function generators reducibles) grob-struct
  90.     (let ((poly (make-epolynomial ring compare-function relation)))
  91.       (push (poly-form poly) reducibles)
  92.       poly)))
  93.  
  94. (defmethod relations ((grob-struct grobner-basis))
  95.   (with-slots (generators reducibles compare-function ring) grob-struct
  96.     (append
  97.      (loop for g in generators
  98.        collect (make-instance 'epolynomial
  99.              :domain ring
  100.              :compare-function compare-function
  101.              :form g))
  102.         (loop for g in reducibles
  103.        collect (make-instance 'epolynomial
  104.              :domain ring
  105.              :compare-function compare-function
  106.              :form g)))))
  107.  
  108. (defmethod reset-grobner-basis ((grob-struct grobner-basis))
  109.   (with-slots (generators undones possibles reducibles) grob-struct
  110.     (setq generators nil undones nil
  111.       possibles nil reducibles nil)))
  112.  
  113. #+Ignore
  114. (defun terms-s-poly (compare-function terms1 terms2)
  115.   (let ((m (max (le terms1) (le terms2))))
  116.     (gterms-difference compare-function
  117.      (gterms-mon-times terms1 (- m (le terms1)) (lc terms2))
  118.      (gterms-mon-times terms2 (- m (le terms2)) (lc terms1)))))
  119.  
  120. ;; The following saves a bunch of consing, but not as much as I would expect
  121. (defun terms-s-poly (compare-function terms1 terms2)
  122.   (let ((m (max (le terms1) (le terms2)))
  123.     (ans-terms (list nil))
  124.     (terms nil)
  125.     sum x y xe ye xc yc new-xe new-ye)
  126.     (macrolet
  127.     ((collect-term (.e. .c.)
  128.        `(progn (setf (rest terms) (make-terms , .e. , .c.))
  129.          (setf terms (rest terms)))))
  130.       (setq terms ans-terms)
  131.       (setq x (red terms1) y (red terms2))
  132.       (setq xe (- m (le terms1))
  133.         xc (lc terms2)
  134.         ye (- m (le terms2))
  135.         yc (- (lc terms1)))
  136.       (loop
  137.        (cond ((terms0? x)
  138.           (cond ((terms0? y) (return (rest ans-terms)))
  139.             (t (collect-term (+ ye (le y)) (* yc (lc y)))
  140.                (setq y (red y)))))
  141.          ((or (terms0? y)
  142.           (%funcall compare-function
  143.                 (setq new-xe (+ xe (le x)))
  144.                 (setq new-ye (+ ye (le y)))))
  145.           (collect-term (+ xe (le x)) (* xc (lc x)))
  146.           (setq x (red x)))
  147.          ((%funcall compare-function new-ye new-xe)
  148.           (collect-term new-ye (* yc (lc y)))
  149.           (setq y (red y)))
  150.          (t (setq sum (+ (* xc (lc x)) (* yc (lc y))))
  151.         (unless (0? sum)
  152.           (collect-term new-xe sum))
  153.         (setq x (red x) y (red y))))))))
  154.  
  155. (defmethod reduce-basis ((grob-struct grobner-basis))
  156.   (with-grobner-operations grob-struct
  157.     (flet ((criterion1 (degree f1 f2)
  158.          (loop for p in generators do
  159.            (when (and (not (eql p f1))
  160.               (not (eql p f2))
  161.               (dominates degree (le p)))
  162.          (unless (member nil undones
  163.                  :test (lambda (x prod)
  164.                      (declare (ignore x))
  165.                      (let ((b1 (second prod))
  166.                            (b2 (third prod)))
  167.                        (or (and (eql f1 b1) (eql p b2))
  168.                            (and (eql f1 b2) (eql p b1))
  169.                            (and (eql p b1) (eql f2 b2))
  170.                            (and (eql p b2) (eql f2 b1))))))
  171.            (return-from criterion1 t))))))
  172.       (let (temp f1 f2 h)
  173.     (reduce-all grob-struct)
  174.     (new-basis grob-struct)
  175.           (loop while undones do
  176.           (setq temp (pop undones))
  177.       (setq f1 (second temp) f2 (third temp))
  178.       (when (and (null (criterion1 (first temp) f1 f2))
  179.              (not (disjoint (le f1) (le f2))))
  180.         (setq h (terms-reduce compare-function
  181.                   (gterms-prim*
  182.                    (terms-s-poly compare-function f1 f2))
  183.                   generators))
  184.         (when (not (terms0? h))
  185.           (setq reducibles nil)
  186.           (setq possibles (list h))
  187.           (setq generators
  188.             (loop for g in generators
  189.               when (dominates (le g) (le h))
  190.                 do (push g reducibles)
  191.               else collect g))
  192.           (setq undones
  193.             (loop for undone in undones
  194.               unless (or (member (second undone) reducibles)
  195.                      (member (third undone) reducibles))
  196.                 collect undone)) 
  197.           (reduce-all grob-struct)
  198.           (new-basis grob-struct)))))))
  199.   grob-struct)
  200.  
  201. ;; This makes sure that all of the polynomials in gneerators and
  202. ;; possibles are AUTOREDUCED.
  203. (defmethod reduce-all ((grob-struct grobner-basis))
  204.   (with-grobner-operations grob-struct
  205.     (let (h g0)
  206.       (loop while (not (null reducibles)) do
  207.         (setq h (terms-reduce compare-function
  208.                   (pop reducibles)
  209.                   (append generators possibles)))
  210.     (unless (terms0? h)
  211.       (setq generators (loop for elt in generators 
  212.                  when (dominates (le elt) (le h))
  213.                    do (push elt reducibles)
  214.                       (push elt g0)
  215.                  else collect elt))
  216.       (setq possibles (loop for elt in possibles
  217.                 when (dominates (le elt) (le h))
  218.                   do (push elt reducibles)
  219.                 else collect elt))
  220.       (setq undones (loop for (nil f1 f2) in undones
  221.                   when (and (not (member f1 g0))
  222.                     (not (member f2 g0)))
  223.                 collect (list (max (le f1) (le f2)) f1 f2)))
  224.       (push h possibles))))))
  225.  
  226. (defmethod new-basis ((grob-struct grobner-basis))
  227.   (with-grobner-operations grob-struct
  228.     (flet ((add-undone (f g)
  229.          (when (e> (le f) (le g))
  230.            (rotatef f g))
  231.          (loop for (nil ff gg) in undones
  232.            do (when (and (eql ff f) (eq gg g))
  233.             (return t))
  234.            finally (push (list (max (le f) (le g)) f g) undones))))
  235.       (setq generators (append generators possibles))
  236.       (loop for g in generators do
  237.     (loop for elt in possibles do
  238.       (when (not (eql elt g))
  239.         (add-undone elt g))))
  240.       (setq possibles nil)
  241.       (setq undones (sort undones (lambda (a b) (e< (first a) (first b)))))
  242.       #+ignore
  243.       (setq generators
  244.         (loop for g in generators
  245.           for h = (terms-reduce compare-function g (remove g generators))
  246.           when (not (terms0? h))
  247.             collect h)))))
  248.  
  249. ;; Reduce terms modulo the current basis
  250. (defun terms-reduce (compare-function terms basis)
  251.   #+ignore
  252.   (format t "~&~%Poly = ~S~%Basis: "
  253.         (le terms))
  254.   #+ignore
  255.   (princ (mapcar (lambda (f) (le f)) basis))
  256.   (let ((again t))
  257.     (loop while again do
  258.       (when (terms0? terms)
  259.     (return nil))
  260.       #+ignore
  261.       (format t "~&Terms = ~S"
  262.         (make-instance 'epolynomial
  263.           :domain (slot-value grob-struct 'ring)
  264.           :compare-function compare-function
  265.           :form terms))
  266.       (loop for g in basis
  267.         do (when (dominates (le terms) (le g))
  268.          (setq terms (gterms-prim*
  269.                   (terms-s-poly compare-function terms g)))
  270.          (return t))
  271.           finally (setq again nil))))
  272.     #+ignore
  273.     (format t "~&Result = ~S~%" (le terms))
  274.     terms)
  275.  
  276. ;; Make poly primitive.  
  277. ;; This isn't really well defined since coefs are in a field.  Idea is
  278. ;; to make the coefficients smaller.  Its really worth avoiding
  279. ;; dividing out a content of 1!!!
  280. #+ignore  ;; Use for integral domains
  281. (defun gterms-prim* (poly) 
  282.   (unless (terms0? poly)
  283.     (let ((coef-domain (domain-of (lc poly)))
  284.       (num-gcd (numerator (lc poly)))
  285.       (den-gcd (denominator (lc poly)))
  286.       1/content)
  287.       ;; Should really use a probabilistic algorithm content algorithm
  288.       ;; here
  289.       (map-over-each-term (red poly) (nil c)
  290.     (if (1? num-gcd)
  291.         (if (1? den-gcd) (return t)
  292.         (setq den-gcd (gcd den-gcd (denominator c))))
  293.         (if (1? den-gcd)
  294.         (setq num-gcd (gcd num-gcd (numerator c)))
  295.         (setq num-gcd (gcd num-gcd (numerator c))
  296.               den-gcd (gcd den-gcd (denominator c))))))
  297.       (unless (and (1? num-gcd) (1? den-gcd))
  298.     (setq 1/content (make-quotient-element coef-domain den-gcd num-gcd))
  299.     (map-over-each-term poly (e c)
  300.       (update-term e (* c 1/content))))))
  301.   poly)
  302.  
  303. ;; Use for fields
  304. (defun gterms-prim* (poly)
  305.   (unless (terms0? poly)
  306.     (let ((1/lc (/ (lc poly))))
  307.       (unless (1? 1/lc)
  308.     (map-over-each-term poly (e c)
  309.       (update-term e (* c 1/lc))))))
  310.   poly)
  311.